## Also need to transform to factor for## setting up the MRF smooth.Irrig_Rev_rice_wheat$District <-as.factor(Irrig_Rev_rice_wheat$a_q103_district)## Now note that not all regions are observed,## therefore we need to remove those regions## from the penalty matrixrn <-rownames(K)lv <-levels(Irrig_Rev_rice_wheat$District)i <- rn %in% lvK <- K[i, i]set.seed(321)library(bamlss)
2 Non-structural: Bivariate gaussian
2.1 Yield analysis
Code
# Bivariate gaussian ## Formulasbiv_gauss_yield<-list( b_grain_yield_ton_per_ha_rice ~1, l_ton_per_hectare ~1, sigma1 ~1, sigma2 ~1, rho ~1)## Model fittingbiv_gauss_yield_model <-bamlss(biv_gauss_yield, family = bamlss:::bivnorm_bamlss, data = Irrig_Rev_rice_wheat,burnin =2000, thin =10, n.iter =12000, sampler =TRUE, nu =1)
## Predictionbiv_gauss_rev_non_struct_pred <-predict(biv_gauss_rev_non_struct, type ="parameter")biv_gauss_rev_non_struct_pred=as.data.frame(biv_gauss_rev_non_struct_pred)
## Predictionbiv_gauss_rev_struct_pred <-predict(biv_gauss_rev_struct, type ="parameter")biv_gauss_rev_struct_pred=as.data.frame(biv_gauss_rev_struct_pred)
4 Modified Cholesky Decomposition Analyses
4.1 Yield analysis
4.1.1 Trial:Non-structural without controls
Code
# Traditional choleskyf_rice_wheat_yield_MRF_corr_endo_yld_non_struct_nocntrl<-list(b_grain_yield_ton_per_ha_rice ~1 , l_ton_per_hectare ~1, lamdiag1 ~1, lamdiag2 ~1, lambda12 ~1)multivariate_geo_sow_MRF_corr_endo_yld_non_struct_nocntrl <-bamlss(f_rice_wheat_yield_MRF_corr_endo_yld_non_struct_nocntrl, family =mvnchol_bamlss(k =2), data = Irrig_Rev_rice_wheat)
## More results -------------------------------nd <-data.frame("District"=unique(Irrig_Rev_rice_wheat$District))# Focusing on the cross-equation correlations## Predict for the structured spatial effects.p_str_multivariate_geo_sow_MRF_corr_endo_rice_y <-predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term ="s(District,id='mrf1')", intercept =FALSE)p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt <-as.data.frame(p_str_multivariate_geo_sow_MRF_corr_endo_rice_y)write.csv(p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt, "tables/p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt2.csv")## And the unstructured spatial effect.p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y <-predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term ="s(District,id='re2')", intercept =FALSE)# Rice sowing spatial equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main =expression(mu(rice_sowing)), title ="Structured spatial effect")
Code
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main =expression(mu(rice_sowing)), title ="Unstructured spatial effect")
Code
# Rice yield spatial equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main =expression(mu(rice_yield)), title ="Structured spatial effect")
Code
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main =expression(mu(rice_yield)), title ="Unstructured spatial effect")
Code
# Wheat sowing equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main =expression(mu(wheat_sowing)), title ="Structured spatial effect")
Code
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main =expression(mu(wheat_sowing)), title ="Unstructured spatial effect")
Code
# Wheat yield spatial equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main =expression(mu(wheat_yield)), title ="Structured spatial effect")
Code
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main =expression(mu(wheat_yield)), title ="Unstructured spatial effect")
Code
# Focusing on the cross-equation correlations## MRF smooth effect.plotmap(India_aoi_sp_bnd,x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District,main =expression(lambda(rice, wheat)), title ="Structured spatial effect")
Code
## Random effects.plotmap(India_aoi_sp_bnd,x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District, main =expression(lambda(rice, wheat)), title ="Unstructured spatial effect")
Code
# Rice sowing equation : Non-linear relationships# s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15)plot(multivariate_geo_sow_MRF_corr_endo, model ="mu1", term ="s(gw_2017)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu1", term ="s(onset_2017)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu1", term ="s(monsoon_onset_dev)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu1", term ="s(median_onset_82_15)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu1", term ="s(sd_onset_82_15)")
Code
# Rice yield equation# s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice)plot(multivariate_geo_sow_MRF_corr_endo, model ="mu2", term ="s(Res_rice_sow)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu2", term ="s(g_q5305_irrig_times_rice)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu2", term ="s(nperha_rice)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu2", term ="s(p2o5perha_rice)")
Code
# Wheat sowing equation# s(harvest_day_rice) + s(gw_2018)plot(multivariate_geo_sow_MRF_corr_endo, model ="mu3", term ="s(harvest_day_rice)")
Code
# (multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(gw_2018)")# Wheat yield equationplot(multivariate_geo_sow_MRF_corr_endo, model ="mu4", term ="s(Res_wheat_sow)")
Code
# Fitted valuesmultivariate_geo_sow_MRF_corr_endo_fitted_values <- multivariate_geo_sow_MRF_corr_endo$fittedmultivariate_geo_sow_MRF_corr_endo_fitted_values <-as.data.frame(multivariate_geo_sow_MRF_corr_endo_fitted_values)# Merge the fitted results to the data and exportIrrig_Rev_rice_wheat_Mult_Res <-cbind(Irrig_Rev_rice_wheat, multivariate_geo_sow_MRF_corr_endo_fitted_values)write.csv(Irrig_Rev_rice_wheat_Mult_Res, "tables/Irrig_Rev_rice_wheat_Mult_Res.csv")summary(lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res))
Call:
lm(formula = sowdate_fmt_rice_day ~ mu1, data = Irrig_Rev_rice_wheat_Mult_Res)
Residuals:
Min 1Q Median 3Q Max
-44.574 -5.397 -0.269 4.986 44.985
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03383 3.41385 -0.01 0.992
mu1 1.00018 0.01771 56.48 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.093 on 4525 degrees of freedom
Multiple R-squared: 0.4135, Adjusted R-squared: 0.4134
F-statistic: 3190 on 1 and 4525 DF, p-value: < 2.2e-16
# Maximum autocorrelationplot(multivariate_geo_sow_MRF_corr_endo_rev, which ="max-acf" , spar =FALSE, lag =200)dev.off()
null device
1
Code
# Traceplots#par(mar = c(4, 4, 4, 1))#plot(multivariate_geo_sow_MRF_corr_endo_rev, which = "samples")#model = "mu1", term = "(Intercept)"# 95% Credible Interval for Predictionsp <-predict(multivariate_geo_sow_MRF_corr_endo_rev, type ="parameter", FUN = c95)p=as.data.frame(p)## More results -------------------------------nd <-data.frame("District"=unique(Irrig_Rev_rice_wheat$District))# Focusing on the cross-equation correlations## Predict for the structured spatial effects.p_str_multivariate_geo_sow_MRF_corr_endo_rice_y_rev <-predict(multivariate_geo_sow_MRF_corr_endo_rev, newdata = nd, term ="s(District,id='mrf1')", intercept =FALSE)p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt_rev <-as.data.frame(p_str_multivariate_geo_sow_MRF_corr_endo_rice_y_rev)#write.csv(p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt_rev, "tables/p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt2_rev.csv")## And the unstructured spatial effect.p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y_rev <-predict(multivariate_geo_sow_MRF_corr_endo_rev, newdata = nd, term ="s(District,id='re2')", intercept =FALSE)# Rice sowing spatial equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y_rev$mu1, id = nd$District, main =expression(mu(rice_sowing)), title ="Structured spatial effect")plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y_rev$mu1, id = nd$District, main =expression(mu(rice_sowing)), title ="Unstructured spatial effect")# Rice yield spatial equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y_rev$mu2, id = nd$District, main =expression(mu(rice_revenues)), title ="Structured spatial effect")plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y_rev$mu2, id = nd$District, main =expression(mu(rice_revenues)), title ="Unstructured spatial effect")# Wheat sowing equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y_rev$mu3, id = nd$District, main =expression(mu(wheat_sowing)), title ="Structured spatial effect")plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y_rev$mu3, id = nd$District, main =expression(mu(wheat_sowing)), title ="Unstructured spatial effect")# Wheat yield spatial equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y_rev$mu4, id = nd$District, main =expression(mu(wheat_revenues)), title ="Structured spatial effect")plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y_rev$mu4, id = nd$District, main =expression(mu(wheat_revenues)), title ="Unstructured spatial effect")# Focusing on the cross-equation correlations## MRF smooth effect.plotmap(India_aoi_sp_bnd,x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y_rev$lambda24, id = nd$District,main =expression(lambda(rice, wheat)), title ="Structured spatial effect")## Random effects.plotmap(India_aoi_sp_bnd,x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y_rev$lambda24, id = nd$District, main =expression(lambda(rice, wheat)), title ="Unstructured spatial effect")# Rice sowing equation : Non-linear relationships# s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15)plot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu1", term ="s(gw_2017)")plot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu1", term ="s(onset_2017)")plot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu1", term ="s(monsoon_onset_dev)")plot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu1", term ="s(median_onset_82_15)")plot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu1", term ="s(sd_onset_82_15)")# Rice yield equation# s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice)plot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu2", term ="s(Res_rice_sow)")plot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu2", term ="s(g_q5305_irrig_times_rice)")plot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu2", term ="s(nperha_rice)")plot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu2", term ="s(p2o5perha_rice)")# Wheat sowing equation# s(harvest_day_rice) + s(gw_2018)plot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu3", term ="s(harvest_day_rice)")# (multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(gw_2018)")# Wheat yield equationplot(multivariate_geo_sow_MRF_corr_endo_rev, model ="mu4", term ="s(Res_wheat_sow)")# Fitted valuesmultivariate_geo_sow_MRF_corr_endo_fitted_values_rev <- multivariate_geo_sow_MRF_corr_endo_rev$fittedmultivariate_geo_sow_MRF_corr_endo_fitted_values_rev <-as.data.frame(multivariate_geo_sow_MRF_corr_endo_fitted_values_rev)# Merge the fitted results to the data and exportIrrig_Rev_rice_wheat_Mult_Res_rev <-cbind(Irrig_Rev_rice_wheat, multivariate_geo_sow_MRF_corr_endo_fitted_values_rev)write.csv(Irrig_Rev_rice_wheat_Mult_Res_rev, "tables/Irrig_Rev_rice_wheat_Mult_Res_rev.csv")summary(lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res_rev))
Call:
lm(formula = sowdate_fmt_rice_day ~ mu1, data = Irrig_Rev_rice_wheat_Mult_Res_rev)
Residuals:
Min 1Q Median 3Q Max
-44.437 -5.403 -0.199 5.055 45.186
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1847 3.4311 -0.054 0.957
mu1 1.0010 0.0178 56.240 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.109 on 4525 degrees of freedom
Multiple R-squared: 0.4114, Adjusted R-squared: 0.4113
F-statistic: 3163 on 1 and 4525 DF, p-value: < 2.2e-16